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SUMMARY 


A developmental  study  was  performed  to  aid  the  AFWL  in  the  iden- 
tification of  soil  properties  from  CIST  test  data.  A mathematical 
algorithm  was  developed  which,  using  prior  estimates  of  the  properties 
and  velocity-time  history  data  from  the  tests,  provided  improved 
estimates  of  the  parameters  in  the  soil  model.  The  algorithm  was 
tested  using  a computational  experiment,  i.e.,  a finite-difference 
code  was  used  to  generate  a set  of  velocity-time  history  responses 
based  on  a predetermined  set  of  parameters.  The  algorithm  was  then 
used  to  determine  a set  of  parameters  based  on  the  data,  and  compari- 
sons were  made  with  the  exact  solution.  After  a number  of  cycles 
through  the  algorithm,  a set  of  parameters  was  derived  which  pro- 
vided satisfactory  matching  of  the  velocity- time  histories. 
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SECTION  I 
INTRODUCTION 

For  several  years,  AFWL  has  been  conducting  a series  of  CIST  (Cylin- 
drical In  Situ  Tests)'*'  tests  to  acquire  data  on  the  shock  response 
of  soils  and  develop  models  which  may  be  used  in  making  predictions 
of  response  of  facilities  to  ground  motion  from  nuclear  weapons 
effects.  Because  of  the  very  high  levels  of  the  shock  in  these  tests, 
simple  linear  soil  models  are  not  adequate.  Consequently,  in  parallel 
to  the  physical  tests,  considerable  effort  has  been  spent  in  support 
of  the  modeling  of  the  tests  using  sophisticated  finite  difference 
codes  such  as  AFTON. 

If  the  assumption  is  made  that  the  mathematical  model  in  the  AFTON 
program  is  valid  then,  with  the  right  values  of  the  soil  parameters 
and  the  proper  representation  of  the  initial  shock,  it  should  be 
possible  to  duplicate  the  test  response  using  the  computer  simulation. 
Unfortunately,  because  of  the  number  of  soil  parameters  involved  and 
because  of  their  interrelationships,  it  has  been  very  difficult  to 
reproduce  the  response  measured  in  tests  by  means  of  the  analyst's 
soil  parameter  estimates  and  computer  simulation.  The  adjustment  of 
the  parameters  to  obtain  the  proper  response  has  involved  a cut  and 
try  method  which  required  a large  effort  by  individuals  who  had  exten- 
sive experience  and  long  exposure  to  prior  data  and  response,  and  has 
resulted  in  very  limited  success. 

The  problem,  as  described,  made  it  very  difficult  to  process  the  large 

quantities  of  CIST  data  and  an  alternative  approach  had  to  be  taken. 

It  was  at  this  point  that  the  J.H.  Wiggins  Company,  which  had  devel- 

2 

oped  a parameter  estimation  algorithm  for  NASA  , proposed  a soil 
parameter  identification  algorithm  which  could  be  used  to  systematize 
and  automate  the  parameter  identification  process  which  had  been  done 
by  hand  up  to  that  date. 

1.  Davis,  S.E.,  Experimental  Data  from  the  Middle  Gust  and  Mixed  Com- 
pany CIST  Events,  Santa  Barbara,  CA,  DNA3151P2,  1 May  1973. 

2.  Collins,  J.D.,  Hart,  G.C.,  Hasselman,  T.K.,  Kennedy,  B. , Statisti- 
cal Identification  of  Structures,  AI.AA  Journal,  Vol.  12,  No.  2, 
pp.  185-190,  February  1974. 
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Since  the  parameter  identification  method  was  untried  in  this  area  of 
mechanics,  it  was  decided  that  a test  problem  be  developed  which  could 
be  used  to  evaluate  the  validity  of  the  identification  methodology.  A 
problem  was  designed  where  the  AFTON  program  was  used  with  a theoret- 
ical model  to  produce  a velocity  time  history  of  response  that  the 
identification  algorithm  was  supposed  to  match.  By  using  the  AFTON 
output  rather  than  the  output  from  a CIST  test,  the  problem  of  model 
validity  and  test  data  validity  did  not  have  to  be  considered.  The 
only  question  was  whether  the  identification  algorithm  could  reproduce 
the  parameters  that  had  been  originally  used  to  obtain  the  velocity 
time  histories. 

The  following  sections  of  this  report  contain  the  results  of  this  ex- 
periment, describe  the  system  identification  algorithm  concept  and 
show  how  the  methodology  is  implemented  on  the  computer. 

The  identification  algorithm,  which  is  referred  to  as  ESP,  is  not  the 
standard  least  squares  type  of  algorithm  which  is  most  frequently  used 
in  parameter  identification.  This  method,  which  will  be  explained 
later  as  a Bayesian  method,  incorporates  the  judgment  of  the  analyst 
into  the  analysis  in  order  to  improve  the  convergence  to  the  proper 
values . 
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SECTION  II 

DEVELOPMENT  OF  THE  ESP  MODEL 
PARAMETER  ESTIMATION 

There  is  a common  problem  that  faces  the  analyst  and  the  experimen- 
talist in  a number  of  fields.  The  analyst,  from  the  physics  of  the 
problem  and  past  experience,  develops  a mathematical  model  to  pre- 
dict the  response  of  the  structure  or  other  media.  Presumably,  given 
a set  of  conditions,  the  model  will  reflect  real  life  behavior.  Un- 
fortunately, when  the  modeling  problem  is  difficult,  the  behavior  as 
measured  in  a physical  test  frequently  does  not  match  satisfactorily 
the  behavior  as  predicted  by  the  model.  This  is  the  dilemma. 

There  are  two  problems  that  can  arise  in  the  modeling  of  the  phenom- 
enon: the  mathematical  model,  and  the  parameters  which  are  used  in  the 
model.  Both  the  model  and  the  parameters  can  be  in  error;  however, 
if  the  model  is  valid  then  it  is  the  simpler  problem  to  find  the  cor- 
rect values  of  the  parameters.  If  the  model  is  invalid,  it  is  very 
unlikely  that  a set  of  parameters  will  be  found  that  will  produce 
satisfactory  results.  In  the  work  performed  in  this  study,  the  assump- 
tion was  made  that  the  model  was  valid.  (This  in  fact  was  accom- 
plished by  developing  pseudo  test  data  from  a mathematical  model  and 
then  determining  whether  the  estimation  method  would  produce  the 
proper  parameters.)  However,  in  most  situations  the  adequacy  of  the 
model  can  still  be  challenged.  The  best  solution  to  this  problem  is 
to  perform  the  parameter  estimation  with  several  models  in  order  to 
find  that  model  and  parameter  set  which  most  adequately  matches  the 
real  situation  (test  data) . 

Assuming  now  that  the  model  is  valid,  consider  the  problem  of  estimat- 
ing the  parameters  to  fit  the  model.  The  most  common  parameter  esti- 
mation method  is  the  least-squares  technique.  With  a simple  least 
squares  analysis,  one  assumes  a linear  relationship  between  the 
variables,  say  x and  y.  Using  test  data  gathered  from  measurements 
of  y as  a function  of  x,  one  is  able  to  estimate  the  slope  and  inter- 
cept, i.e.,  the  parameters,  for  the  linear  equation.  With  simple 
linear  models  and  a quantity  of  data  points,  it  is  a relatively  easy 
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task  to  obtain  the  appropriate'  values  of  the  parameters.  However, 
when  the  model  becomes  complex  with  many  parameters  and  a sparsity 
of  data,  the  least-squares  method  becomes  inadequate  for  solving 
the  problem.  Variations  on  the  least-squares  method,  such  as  weighted 
least-squares  and  maximum  likelihood  are  also  limited  in  solving  this 
multiple-parameter,  limited-data  problem.  Consequently,  an  alterna- 
tive approach  had  to  be  found  which  could  incorporate  prior  judgments 
or  knowledge  of  model  parameters  to  aid  the  parameter  identification 
problem.  This  alternative  approach  involves  the  use  of  Bayes’  Rule, 
which  is  explained  in  the  following  section. 

BAYES ’ RULE 

The  fundamental  concept  embodied  in  the  parameter  estimator  used  in 
this  work  is  Bayes1  rule  from  probability  theory.  This  formula, 
originally  a mathematical  curiosity,  has,  during  the  past  twenty 
years,  become  a basis  for  modern  decision  theory. 

The  equation  states: 

P(BlA)P(A) (1) 

P(B|A)P(A)  + P(B|A)P(A) 

where  P(a]b)  is  the  probability  (conditional)  that  event  A will  occur 
given  the  occurrence  of  B,  and  P(A)  = 1 - P (A) . The  equation  can  also 
be  written  in  the  probability  density  form 

p(*ly>  - ,-p(g|%!fu)dt  p(t)  <2) 

The  probability  P (A)  in  equation  (2)  can  be  considered  as  the  first 
estimate  of  the  probability  of  event  A.  P(B|a)  is  the  likelihood  of 
the  particular  B occurring  given  that  A occurred.  The  P(A|B)  is  the 
revision  of  the  probability  of  A given  the  particular  B.  In  words, 
the  equation  is: 

Prior  (first  estimate)  * Likelihood  (comparison  of 
data  and  estimate)  -*■  Posterior  (revised  estimate) 


P (A  j 
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The  significance  of  the  formula  is  that  one  can  guess  a set  of  param- 
eters (the  prior) , and  then  make  a measurement.  The  likelihood  of 
that  particular  measurement,  given  the  prior,  is  then  multiplied  by 
the  prior  and  normalized  to  give  a new  estimate  of  the  set  of  param- 
eters (the  posterior) . 

Bayes'  rule  is  the  only  formal  place  in  probability  and  statistical 
theory  where  human  judgment  can  be  combined  with  measurement  to  pro- 
duce a better  answer.  This  characteristic  is  particularly  desirable 
when  measured  data  are  meager  and  past  experience  can  be  used  to 
guide  it  in  the  right  direction.  If  there  is  a great  abundance  of 
data,  Bayes'  rule  usually  loses  its  significance. 

In  the  application  to  soil  model  parameters,  the  analyst  is  seeking 
a "posterior  model"  of  the  hydrostat,  the  failure  surface,  and 
Poisson's  ratio.  Hopefully  this  posterior  model  will  reflect  the 
actual  dynamic  behavior  of  the  soil  in  the  test  and  can  be  considered 
to  be  the  best  final  estimate  of  the  actual  soil  model.  To  get  the 
posterior  model,  the  analyst  first  estimates  the  parameters  of  the 
hydrostat  and  the  failure  surface  (the  prior)  and  then  within  the 
estimation  algorithm  establishes  the  probability  of  obtaining  the 
particular  test  results  with  this  model  (the  likelihood) . 

In  this  case,  the  estimated  soil  model  parameters  (hydrostat,  failure 
surface,  and  Poisson's  ratio)  and  their  uncertainties  form  the  prior. 
The  procedure  in  the  following  development  will  show  how  the  prior 
statistical  mode],  is  combined  with  the  test  measurements  to  obtain 
the  posterior  model.  The  posterior  model  is,  ideally,  the  identified 
system. 

EXPANSION  AROUND  A REFERENCE  VALUE 

The  shock  response  of  a soil  system  may  be  represented  like  any  other 

mathematically  well-behaved  functions  in  terms  of  a Taylor's  series. 

For  m mathematical  functions  (f ..  , f-,....f  ) of  n variables  (s.,  a- , 

1 z m i 

...  a ) this  yields  in  vector  form: 
n 1 


_ a_ 
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r 


(f(a)}=<  . 


fl  (al' 

a2 ' 

....  a ) 
n 

f2  ^al ' 

a2 ' 

....  ot  ) 
n 

fm(V 

a2. 

• • . . ct  ) 
n 

(3) 


In  terms  of  soil  mechanics  this  means  that  for  a given  set  of  initial 
conditions  the  state  of  stress,  the  velocity  and  the  strain  of  the 
material  are  functions,  (f(a)},  of  the  constitutive  parameters  of  the 
material,  (a^,  a? , a^,...}.  For  example,  the  peak  particle  veloci- 
ties in  the  soil  at  given  points,  for  a set  of  initial  conditions, 
are  functions  of  quantities  such  as  mass  density,  Poisson's  ratio 
and  wave  speeds.  For  this  example  the  number  of  responses  is  the 
total  number  of  locations  at  which  the  peak  velocity  is  measured. 


The  system  of  equations  can  be  expanded  into  a Taylor's  series 

around  a reference  vector  {a  } 

P 


or 


1/P 

‘2,p 


{ f («)  } = { f (oip)  } + £■ 


3f  (a) 
3 (a) 


{ a - a } + 0({  a - a }{a-a  }T)  (4) 

P v P P ' 


where  0 denotes  the  order  of  magnitude  of  the  truncation  error 


{ f (a  ) } 


= < 


a 2 , p ' 

...  a ) 

n , p 

£2<“l,p' 

a2 , p ' 

...  a , 

n,  p) 

a2  , p ' 

. . . a ) 

n , p 
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[Wl  = 

(_  3 (a)J 


3fl 

3 f 1 

3 f 1 

3a^ 

3a2  ' * ' 

3a 

n 

3 f 2 

3 f 2 

3 f 2 

3a^ 

3a2  • • * 

3a 

n 

3 f 

m 

3f 

m 

3 f 

m 

3a^ 

3a2 

8a 

n 

The  Taylor's  series  expansion  shown  here  is  valid  when  f(o.)  possesses 
continuous  first  and  second  derivatives.  Indeed  the  truncation  error 
represented  by  the  term  ©({a-a^}  {a  - a^}7)  is  zero  when  the  sta- 
tistical model  is  linear.  It  is  fortunate  that  most  functions  govern- 
ing physical  phenomena  are  expandable  into  a Taylor's  sexies,  thus 
permitting  the  above  operations  to  be  carried  out. 

It  should  be  noted  that  the  shock  response  of  a soil  system  is  always 
a nonlinear  function  of  the  soil  state  vector,  {a}.  The  assumption 
of  the  truncated  Taylor's  series  shown  above  is  appropriate  when  the 
changes  in  the  soil  parameter  {a-a^}  are  small.  Obviously  as  these 
changes  approach  zero,  the  truncation  error  will  approach  zero  even 
faster  and  the  formulation  given  by  the  truncated  series  becomes,  in 
the  limit,  exact.  Thus  the  repeated  application  of  a seemingly  linear 
statistical  model  can  be  employed  to  solve  a nonlinear  statistical 
model,  by  ignoring  the  truncation  error. 


By  dropping  the  truncation  term,  equation  (4) 
the  approximation 

will  be  simplified  to 

{Af}  = 

(5) 

where 

{ A f 1 = 

i ui  J 

(f (a)  - f (a  ) } 
F 

and 

{Aa}  = 

{ a - a } 

p 
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If  the  number  of  responses,  m,  is  the  same  as  the  number  of  param- 
eters, n,  and  if  the  matrix 


faf  (o)l 
L 8 (a)  J 

is  nonsingular,  equation  (5)  can  be  readily  used  to  solve  for  {Aa}. 
However,  m and  n are  usually  not  equal  and  estimation  techniques 
must  be  resorted  to  in  order  to  determine  {Aa}. 

STATISTICAL  DEFINITIONS 

Before  outlining  the  statistical  identification  process,  it  is  imper- 
ative that  statistical  concepts  be  reviewed. 

E(x)  or  x is  the  expected  or  mean  value  of  the  random  variable  x. 
Since  E is  a linear  operator,  the  expected  value  of  a vector  {x} 
can  be  readily  computed  to  be 


E 


fE(xl 
E (x2 


(6) 


The  variance  (or  standard  deviation  squared)  of  a random  variable 
is 


Var (x± ) 


(7) 


The  covariance  of  two  random  variables  is 


Cov (x . x . ) = a 

i : x . x 


xi>  (xj 


p . .0.0, 

1}  1 3 


(S) 


where  p . 
variables 


is  defined  as  the 

x.  and  x.  and  can 
i 3 


correlation  coefficient  for  random 
be  shown  to  have  a range  of  -1  <_p^j  <_1. 
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The  variances  and  covariances  of  a vector  are  obtained  by  taking 
the  expected  value  of  the  outer  product  of  the  vector 

E ^{x  - x}{x  - x}T)  = E ^{Ax}{Ax}*  ) = [sxx] 


’E  K - 

^>2) 

E 

k(xL  - Xx)  (x7  - x2)) 

E ((Xl  - 

j1hx2-;2 

») 

e((x2  -x2)2) 

• 

. 

• 

• 

- 

• 

• 

“ 2 

a 

X1 

a 

X1X2  • 

2 

CTX  Pi2°x  ax 

Xi  X1  X2 

-] 

a 

X1X2 

2 

a 

x2 

2 

p,  -,o  a a 

12  X1  x2  x2 

• 

• 

• • 

• • 

— 

If  the  vector  {Az}  = {z  - z}  is  related  to  {Ax}  by  the  linear 
transformation 

{ Az}  = [c] {Ax} 


then 


[SI  = E ({  A z } { A z } T ) 
zz  v 

= E ([c]  {Ax}{Ax}T[c]T) 
= [c] E ({ Ax} { Ax)T){ c}T 

= [=nsXJtnc]T 


(10) 


(ID 
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DERIVATION  OF  STATISTICAL  SYSTEM  IDENTIFICATION  EQUATIONS 

Earlier  in  this  section  an  expression  was  developed  which  related 
the  response  of  a soil  system  to  its  soil  parameters.  The  relation- 
ship will  now  be  used  to  solve  the  inverse  problem  of  determining 
the  values  of  the  soil  parameters  when  the  response  is  specified. 

The  derivation  procedure  used  is  based  on  that  of  a modified 
weighted  least  squares. 

The  cases  when  no  errors  are  present  in  the  experimental  error  will 
be  treated  first.  If  the  mathematical  function  f (a)  is  the  true 
representation  of  the  physical  phenomenon,  then  the  measured  response 
fr  will  be  equal  to  f(a)  appearing  in  equation  (3).  Under  these  con- 
ditions, equation  (5)  can  now  be  recast  as 

{F}  = [T]{A}  (12) 

where 

{F}  = {fr  - f (a  ) } 

and  {A}  = (a-  a } 

P 

It  is  important  to  determine  the  statistical  behavior  of  the  random 
variables  in  equation  (12) . It  can  be  assumed  that  the  prior  assumed 
values  of  the  soil  parameters  are  distributed  with  mean 

E ( { a })  = {a}  (13) 

P 

The  covariance  matrix  of  {a}  which  must  be  specified  by  the  soil- 
analyst  is 


[S  ] = E ({a  - a } { a - a }T) 
a a v o d 


(14) 


From  equation  (13),  the  mean  and  covariance  matrix  of  A-  are 


and 


E ( { A } ) = {0} 


[Saa]  = E((AHA}a)  - [Saa] 


(15) 

(16) 
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If  in  equation  (12) , (A)  is  defined  by  a multivariate  normal  distri- 
bution, then  {F}  through  the  linear  transformation  is  also  normally 
distributed  with  mean 

E ( { F } ) = 0 


and  covariance 


[Spp]  = E({F){F}T)  = E(  [T]  {aHa}T{T}T) 


" [T] [Saa]  tT] 


(17) 


The  above  relation  must  now  be  amended  to  account  for  the  error  in 
the  measurements,  {e}.  These  errors  are  independent  of  the  uncertainty 
resulting  from  the  soil  parameter  estimates,  but  do  not  appear  in 
equation  (12).  To  account  for  this  measurement  error,  the  vector  {e} 
is  added  to  the  right-hand  side  of  equation  (12) , giving 


iF)  = f T ] { A } + {e} 


(18) 


The  vector  ( £ } being  unbiased  has  a zero  mean  and  a covariance  (Sc  . ] 
specified  from  the  test  measurements.  The  mean  of  {F} is  still  zero, 
but  the  covariance  is 


Isff^ 


The  covariance  of 
ment  is 


E({F} {F}T) 

e(([t]{a}  + {£})([t](a1  + { e } ) T ) 

IT] [Saa] [T]T  + [SeE]  (19) 

^F}  and  (A)  which  is  required  later  in  the  develop- 


[SFA1  = E ( { F } { A *) 

= E (( CTI f A}  + ( :}) (A  T ) 

= E([T|(AKAlT  + {e}{A/T) 
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■ tTHS^l  = [T][Saa]  (20) 

m 

where  the  expected  value  of  {c}{a}  is  zero  since  the  vectors  {e}  and 
{A}  are  statistically  independent. 

The  purpose  of  this  process  is  to  find  a best  estimator  of  {A}  based 
on  the  measured  responses  { f ^ } and  the  prior  estimates  of  the  soil 
parameters  { ot^ } along  with  their  associated  errors.  We  are  therefore 
seeking  an  equation  of  the  form 

{A*}  = [G] {F}  (21) 

The  matrix  [G]  will  now  be  defined  in  such  a way  as  to  minimize  the 
variance  of  the  difference  between  the  true  value  of  {A*}  and  its 
estimated  value.  Quantitatively,  the  equation 

[Q]  = E [ { A*  - A} {A*  - A}T]  (22) 


is  minimized  with  respect  to  [G] . First  substitute  equation  (21) 
into  (22)  yielding 


[Q] 


e[([G]{F}  - { A } ) ( [ G ] { F } - {A})T] 
[G][Sff][G]T  - [G][SFA]  - [SFA]T[G]T  + [Saa] 


(23) 


Next,  taking  the  variation  of  this  with  respect  to  [G] , we  obtain 
0=  [5G]  ([SFF]  [G]T  - [SFA])  + ( [G]  [SFp]  - [Sfa]T)[5G]T 


(24) 


This  yields 


[G]  - [Sfa]T|Sff]-1 

(a*)  = [sFAlT[sFFr1{Fi  (25) 

which  after  substitution  from  equations  ( 12 ) , (19)  and  (20)  becomes 
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{a*}  = {ap}  + [Saa]  tT]T([T]  [Saa]  [T]T+  [S£c])_1  {fr-f(<*  )i 


The  covariance  of  the  estimate  of  {a*}  is 


(26) 


[S*  ] 

aa 


= E [ { A*  - A} (A*  - A 


^ SAA^  “ ^SFA^  ^SFF^  ^ SFA^ 


= - tSM][T]T([T][Saa][T]T+  [S^])'1  [T][Safl(] 
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ESTIMATOR  OPERATION 

Figure  1 describes  the  sequence  of  operations  used  in  the  soil  model 
parameter  identification  process.  Here  { t } and  {u}  are  the  response 
variables  contained  in  { f } . If  the  relation  between  the  residual  in 
the  observation  vector  { F } and  the  state  vector  {a}  were  indeed  linear, 
the  method  would  yield  a best  estimate  in  one  step  without  iteration. 
However,  the  responses  that  are  computed  are  usually  linked  to  the 
state  vector  in  some  nonlinear  fashion.  Consequently,  a solution  must 
be  obtained  by  repeated  iteration  using  the  successive  approximation 
method. 

As  shown  in  figure  1,  the  experimentalist  supplies  the  ESP  (Estima- 
tion of  Soil  Parameters)  program  with  measured  responses  and  uncer- 
tainties from  CIST  experiments.  The  analyst  supplies  the  program 
with  an  estimate  of  the  state  vector  of  soil  parameters  along  with 
their  uncertainties.  The  procedure  is  then  cycled  between  the  AFTON 
and  ESP  programs  until  convergence  is  obtained. 

By  far  the  most  lengthy  chain  in  the  procedure  just  described  is  the 
calculation  of  the  partial  derivative  matrix.  It  should  be  noted 
that  it  is  the  number  of  parameters  rather  than  responses  that  dic- 
tates the  computational  time  involved  for  this  calculation.  Thus, 
for  example,  the  number  of  computations  used  for  the  12  parameters 
studied  here  was  13  (1  for  the  unperturbed  state  and  12  for  the  in- 
dividually perturbed  states) . 
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Figure  1.  Sequence  of  Operations  of  the 
Soils  Identif ication  Process 
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OBSERVATIONS  ON  THE  ESTIMATOR 

The  procedure  used  in  the  ESP  program  continuously  updates  the  prior 
estimate  while  using  the  uncertainties  of  the  initial  prior  estimate 
as  originally  specified  by  the  analyst.  The  procedure  can  be  im- 
proved by  minimizing  the  weighted  square  of  the  difference  between 
the  initial  prior  value  and  the  best  estimate.  In  such  a sequence 
three  state  vectors  would  be  involved;  the  initial  prior  estimate, 
the  iterative  estimate  and  the  unknown  best  estimate.  The  sequence 
would  iterate  until  the  iterative  estimate  and  the  unknown  best  esti- 
mate are  equal.  Obviously,  the  effect  of  this  suggested  change  to 
the  ESP  program  must  be  studied  before  any  conclusions  can  be  drawn. 

It  should  be  noted  that  the  value  of  the  best  estimate  obtained  in 
the  ESP  program  not  only  depends  on  the  experimental  response  and 
prior  estimates,  but  also  on  their  associated  uncertainties.  When 
the  uncertainty  of  the  experimental  response  is  large,  the  value  of 
the  best  estimate  will  be  close  to  the  prior.  Conversely,  when  the 
uncertainty  of  the  prior  estimate  is  large,  the  value  of  the  best 
estimate  will  be  such  that  the  calculated  response  will  be  close  to 
the  experimental  response.  Thus,  it  is  important  that  the  uncertain- 
ties be  chosen  with  care. 
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SECTION  III 

TEST  OF  THE  ESP  PROGRAM 


BACKGROUND 


In  order  to  obtain  in  situ  material  properties  for  soils  and  rock  for 
use  ' n predicting  crater  and  ground  shock  due  to  nuclear  attack,  AFWL 
has  conducted  CIST  (Cylindrical  In  Situ  Test)  experiments  in  the  U.S., 
Alaska  and  in  the  Pacific.  In  these  tests  PETN  explosive  (more 
commonly  known  as  Primacord) , with  a weight  of  5 pounds  per  linear 
foot,  is  inserted  into  a 2-foot  diameter  hole  to  a depth  ranging 
from  40  to  80  feet.  The  soil  layers  around  the  hole  are  instrumented 
with  accelerometers  which  measure  acceleration  as  a function  of  time 
in  the  CIST  event.  A typical  example  of  a velocity  waveform  gener- 
ated by  integrating  the  acceleration  with  time  is  depicted  in  figure 
2.  Velocity  response  is  of  particular  interest  because  it  can  be 
readily  interpreted  and  compared  to  the  predicted  waveform  generated 
by  the  mathematical  finite  difference  model.  For  a given  location, 
the  significant  responses  which  can  be  used  to  describe  the  wave- 
form are  the  arrival  rime,  bow  time,  bow  velocity,  peak  time,  peak 
velocity,  and  time  to  1/2  peak.  The  arrival  time  is  arbitrarily  de- 
fined to  be  the  time  when  the  velocity  rises  to  5%  of  the  peak 
velocity.  (The  random  fluctuations  shown  prior  to  arrival  represent 
background  noise.)  The  bow  point  is  defined  as  the  maximum  perpen- 
dicular distance  of  the  waveform  from  a line  connecting  the  arrival 
point  to  the  peak  point.  The  1/2  peak  time  is  defined  at  that  time 
when  the  velocity  falls  to  half  the  peak  velocity. 

The  simulation  of  ground  response  at  AFWL  is  executed  by  means  of  a 
set  of  computer  codes  known  as  AFTON.  This  name  specifies  programs 
used  to  solve  transient  continuum  motion  problems.  To  date,  three 
AFTON  codes  have  been  developed:  AFTON  1 is  used  to  solve  one-dimen- 

sional problems,  AFTON  2P  simulates  flow  in  plane-symmetric  systems, 
and  AFTON  2A  is  used  to  solve  flow  in  axisymmetric  systems  where 
radius  and  axial  position  are  the  space  variables.  AFTON  2A  is  the 
program  which  models  the  CIST  experiment  described  above. 
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PEAK  VELOCITY 


Figure  2 . Representative  Waveform 
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Briefly,  AFTON  2A  expresses  the  conservation  of  mass,  momentum  and 
energy  in  a set  of  explicit  finite  d_fference  approximations  that 
march  forward  in  time.  Due  to  considerations  of  stability  and  con- 
vergence, approximations  of  this  form  have  an  upper  limit  to  the 
size  of  their  time  increment. 

In  this  research  program,  twelve  soil  parameters  were  stipulated 
by  AFWL  as  uniquely  describing  the  modeling  of  a particular  soil. 
These  parameters  are  associated  with  the  hydrostat  and  failure 
surface  of  the  soil. 

The  hydrostat  (figure  3)  specifies  the  pressure  as  a function  of 
excess  compression , (p/pa)-l , where  p denotes  density  and  pQ  denotes 
initial  density.  The  slope  of  the  pressure-excess  compression  curve 
is  the  bulk  modulus,  BK,  and  is  calculated  by  the  relation 


where  v and  c are  Poisson's  ratio  and  the  confined  sound  speed, 

respectively.  AFWL  approximates  the  hydrostat  by  assuming  four 

regions  of  behavior.  The  loading  and  unloading  curves  for  the  first 

three  regions  are  assumed  to  be  piecewise  linear.  The  first  region, 

where  0 £ y £ y^,  is  known  as  the  seismic  toe  and  is  characterized 

by  disturbances  so  small  that  both  the  loading  and  unloading  curves 

essentially  follow  the  same  line  with  v = v and  c = c . In  the 

u u 

second  region  where  y^  <_  y £ progressive  fracturing  of  the  soil 

induces  softening  and  produces  a smaller  bulk  modulus  with  v = vL  and 

c = cT  during  loading.  At  the  same  time,  the  unloading  curve  has  a 
"l 

slope,  BK,  identical  to  that  in  the  first  region.  For  the  third 

region  where  ^2  V £ t,,  the  soil  regains  its  stiffness  or  softens 

further  and  the  bulk  modulus  increases  or  decreases  over  that  of  the 

second  region  with  v = and  c = c^  . (Obviously,  if  3KL2  <_  BKL1 , 

C,  < CT  . ) Again,  as  in  the  first  ?wo  regions,  the  unloading  curves 
L2  L1 

for  the  third  region  follow  the  same  slope.  The  final  region  is  that 
of  lockup  and  is  characterized  by  a more  elastic  behavior  of  the  soil. 
In  this  reaion,  the  bulk  modulus  is  initially  determined  by  v = v and 
c = c , but  the  slope  rises  exponentially  and  eventually  attains  an 
asymptotic  value  of  10  psi.  However,  in  unloading,  the  curve  fellows 
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! 
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the  loading  curve  back  down  to  y = ^ and  then  follows  the  same 
unloading  slope  as  in  the  other  regions.  Therefore,  if  cu  and  c„ 
are  not  equal,  there  will  be  a discontinuity  in  the  bulk  modulus  at 
y = For  this  reason,  c^  and  c z are  usually  set  equal  to  each 

other. 

The  failure  surface  (figure  4)  specifies  the  upper  limit  of  the 
square  root  of  the  second  invariant  of  the  stress  deviator  tensor, 
\/j^  , as  a function  of  the  hydrostatic  pressure  defined  as  one- 
third  of  the  first  invariant  of  the  stress  tensor,  J . Here  AFWL 
divides  the  failure  surface  into  two  regions.  In  the  first  part, 
the  function  rises  linearly  from  Y1  at  a slope  of  SI  to  a value  of 
VM1 , called  the  Von  Mises  limit.  In  the  second  region,  the  function 
is  taken  as  a constant,  VM1. 

The  12  soil  parameters  governing  the  mathematical  model  are  vT , v , 

Li  U 

cT  , cT  , c , c , p. , p_,  y,,  Yl,  SI,  and  VM1.  For  a given  initial 

Li-.  Li~  U Z L A J 

condition  at  the  circumference  of  the  cylindrical  hole,  these  para- 
meters will  determine  the  response  of  the  soil. 

THE  EXPERIMENT 

In  order  to  qualify  the  applicability  of  the  ESP  program  to  CIST 
and  AFTON,  a computational  experiment  was  formulated  by  AFWL  and 
CERF.  In  this  experiment,  a set  of  12  values  was  defined  for  the 
soil  parameters  and  the  soil  response  was  generated  by  AFTON  in 
terms  of  the  aforementioned  six  responses  (two  velocities  and  four 
time  responses)  at  radii  of  3,  5,  and  8 feet  for  a total  of  18 
responses.  For  this  calculation  the  boundary  condition  simulating 
the  explosion  was  a pressure  pulse  with  a time-dependent  exponential 
decay  (p  = 7000  e psi]).  The  18  responses  were  then  supplied 

to  the  J.  H.  Wiggins  Company  as  experimental  measurements,  along 
with  the  estimated  error  in  these  measurements.  Then  an  initial 
estimate  of  values  of  the  soil  model  parameters  was  made  by  AFWL 
(with  no  knowledge  of  the  actual  value  of  the  parameters)  and 
supplied  to  the  Wiggins  Company  as  the  prior  estimate  of  soil  para- 
meters along  with  the  estimated  error  in  these  parameters.  The  task 
of  the  Wiggins  Company  was  to  determine  best  estimates  of  the  soil 
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parameters  that  yielded  responses  that  were  statistically  within  the 
bounds  of  the  measured  responses,  by  using  the  prior  estimate  suppli- 
ed by  AFWL. 

For  this  computational  experiment  the  mathematical  model  was  assumed 
to  be  one-dimensional  with  respect  to  radius.  The  computed  values 
and  the  perturbations  (required  for  the  partial  derivative  approxi- 
mations) were  computed  in  one  AFTON  13-layer  run.  For  this  run  each 
layer  was  so  thick  (14  feet  or  more)  that  the  ground  responses  were 
essentially  one-dimensional. 

In  the  first  phase,  the  statistical  identification  concept  which 

3 

had  been  used  successfully  in  the  NASA  MOUSE  study  was  altered  to 
incorporate  partial  derivative  approximations  as  generated  by  AFTON. 
The  resultant  program  has  been  dubbed  the  ESP  (Estimation  of  Soil 
Parameters)  program. 

The  parameters  of  the  ESP  program  are  shown  in  Table  1.  The  con- 
straints for  these  parameters  are: 


1. 

0 < ci.^  < 0.5 

2. 

0 < <*2  < 0 . 5 

3. 

0 <a3  <a5 

4. 

0 <a4  <a5 

5. 

a6  - a5 

6. 

a5  > 0 

7. 

ct7  > 0 

8. 

a8<a7 

9. 

a9>U2 

10. 

aio  ~ 0 

11. 

au>0 

12. 

a12  -a10 

3.  Collins,  J.  D.,  Hart,  G.  C. , Hasselman,  T.  K. , Kennedy,  B. , 

Model  Optimization  Using  Statistical  Estimation,  J.  H.  Wiggins 
Company,  California,  1973 
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Table  1.  PARAMETERS  OF  THE  ESP  PROGRAM 


Parameter 


Alternate  Name 


'Ll 


'L2 


'1 


*10 


°1 

VM1 


Ml 


12 


Units 

Dimensionless 

Dimensionless 

ft/sec 

ft/ sec 

ft/sec 

ft/sec 

psi 

psi 

Dimensionless 

psi 

Dimensionless 

psi 
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where 


P 


2 


1.5  x 106 
(30.48) 2 


(29) 


These  constraints  arise  either  from  classical  elasticity-plasticity 
or  from  assumptions  in  the  finite  difference  model. 

A schematic  of  the  computational  sequence  is  presented  in  figure  5. 
After  AFWL  defined  the  true  parameters  and  provided  the  Wiggins 
Company  with  the  initial  estimate  of  the  parameters  along  with  the 
correct  response,  the  sequence  entered  a computational  loop.  In 
this  loop,  the  Wiggins  Company,  using  the  statistical  identification 
process,  calculated  a modified  set  of  parameters  and  tested  them  for 
convergence.  If  convergence  had  not  been  obtained,  AFWL  re-ran 
AFTON  to  obtain  the  response  and  partial  derivative  matrix  associ- 
ated with  the  modified  set  of  parameters. 

RESULTS 

The  ESP  program  cycled  the  test  case  through  six  iterations  and  the 
results  indicated  the  pertinent  soil  parameters  to  be  converging  to 
the  statistical  limits  of  the  true  parameter  vector.  Table  2 
summarizes  the  results  obtained  for  the  responses  for  each  iteration. 
This  convergence  is  shown  graphically  in  figure  6,  where  the  re- 
sponses are  depicted  for  stations  at  3,  5 and  8 feet  from  the  point 
of  detonation.  As  the  waveforms  show,  the  comparison  between  results 
for  the  fourth  iteration  and  the  actual  response  is  especially  good 
for  the  frontal  portion  of  the  response.  (Although  a total  of  6 
iterations  were  carried  out,  the  waveforms  for  the  last  2 iterations 
were  received  at  the  Wiggins  Company  after  the  illustrations  had  been 
completed.)  This  is  to  be  expected  since  most  of  the  measurements 
provided  to  the  ESP  program  were  taken  from  this  part  of  the  wave- 
form. 

Table  3 summarizes  the  results  for  the  soil  parameters  for  each 
iteration.  The  last  row  indicates  the  number  of  the  constraints 
violated.  violated  its  constraint  in  the  first  and  third 
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iterations.  In  each  case,  was  adjusted  to  the  limit  of  its 
constraint.  Initially,  this  constraint  was  0.5,  but  this  was 
changed  to  0.43  for  purposes  of  calculating  a physically  tenable 
forward  difference  for  the  partial  derivative  approximation  with 
a perturbation  of  10%.  For  all  iterations,  a_  exceeded  even 

D 0 

after  adjustment.  After  some  consideration,  it  was  determined  that 
this  was  because  the  estimation  technique  was  underestimating  lock- 
up strain  but  compensating  for  it  by  making  a,.  greater  than  a^.  The 
next  constraint  violated  is  that  where  > a g (p^  > p^)  in  the  first 
two  iterations.  However,  the  modified  values  in  iteration  3 and  4 
obey  the  constraint.  Also,  it  is  noted  that  the  values  obtained  for 
ot^  were  near  zero.  This  is  plausible  since  the  actual  yield  surface 
is  a constant  and  the  slope  parameter  should  have  no  influence  upon 
it. 

In  order  to  calculate  the  partial  derivatives,  a forward  difference 
approximation  was  employed  where  the  independent  variable  was  per- 
turbed by  10  per  cent.  Since  it  was  believed  that  such  an  approxi- 
mation may  be  inaccurate  for  nonlinear  problems,  the  perturbation 
was  changed  to  1 per  cent  for  the  fourth  iteration.  However,  it 
was  found  that  the  subsequent  perturbation  was  too  small  for  some 
of  the  time  responses  to  show  a significant  change.  Nonetheless, 
the  resulting  state  vector  of  parameters  for  this  iteration  indicates 
a marked  improvement  over  those  obtained  in  the  third  iteration. 

Figure  7 graphically  compares  the  convergence  of  the  fourth  iter- 
ation toward  the  statistical  bounds  of  the  actual  response.  As  seen, 
the  comparison  is  poorer  as  one  advances  toward  the  negative  phase  of 
the  response. 

After  four  iterations  AFWL  compared  the  waveforms  generated  with  the 
actual  curves.  The  comparison  indicated  that  the  process  was  pro- 
ceeding satisfactorily  toward  convergence  and  the  computational 
experiment  was  deemed  to  be  successfully  completed. 

While  this  report  was  being  completed  two  more  iterations  were 
informally  run  in  conjunction  with  AFWL.  Figure  8 shows  a comparison 
of  the  hydrostat  obtained  from  the  sixth  iteration  with  the  initial 
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and  experimental  hydrostats.  (The  hydrostats  for  the  fourth  and 
sixth  iterations  are  substantially  the  same.)  As  seen,  the  system 
identification  process  is  moving  the  iterative  hydrostat  toward  the 
experiment  curve,  with  good  results  observed  in  the  leading  curve 
and  lock-up  point. 

Figure  9 shows  a comparison  of  the  yield  surface  obtained  from  the 
sixth  iteration  with  the  initial  and  experimental  failure  surfaces. 
As  seen,  the  iterative  procedure  has  moved  the  result  for  the  sixth 
iteration  quite  close  to  the  actual  curve,  while  at  the  same  time 
flattening  it. 
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SECTION  IV 

CONCLUSIONS  AND  OBSERVATIONS 

(1)  The  process  of  statistical  identification  has  been  success- 
fully demonstrated  for  the  estimation  of  soil  model  parameters 
in  material-response  models. 

The  use  of  the  ESP  program  in  conjunction  with  a computational 
experiment  devised  by  AFWL  has  demonstrated  the  capability  of 
the  method  to  yield  an  estimate  of  model  parameters  that  gener- 
ates a response  within  the  uncertainty  of  the  experimental 
response. 

(2)  A better  estimate  can  be  obtained  by  using  more  experimental 
data  for  each  CIST  event. 

A larger  amount  of  experimental  data  will  give  a more  reliable 
best  estimate  and  also  impart  a response  that  closely  approxi- 
mates the  experimental  waveform. 

(3)  Times  should  be  judiciously  used  as  response  variables. 

Care  must  be  exercised  in  including  times  as  response  variables 
because  of  the  relatively  large  discretization  errors  associ- 
ated with  them  in  the  finite  difference  modal.  More  research 
is  required  before  any  final  conclusions  can  be  drawn. 

(4)  The  statistical  identification  process  should  be  tested  on  a 
real  material. 

Since  the  computational  experiment  has  verified  the  application 
of  the  statistical  identification  process  to  material-response 
models,  it  is  suggested  that  the  process  be  tested  on  a real 
material  using  measured  data. 
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